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Abstract. We study the out-of-equilibrium dynamics of the fully-frustrated XY 
model. At equilibrium, this model undergoes two phase transitions at two very 
close temperatures: a Kosterlitz-Thouless topological transition and a second- 
order phase transition between a paramagnetic phase and a low-temperature phase 
where the chiralities of the lattice plaquettes are anti-ferromagnetically ordered. 
We compute by Monte Carlo simulations two-time spin-spin and chirality-chirality 
autocorrelation and response functions. From the dynamics of the spin waves in 
the low temperature phase, we extract the temperature-dependent exponent r\. 
We provide evidences for logarithmic corrections above the Kosterlitz-Thouless 
temperature and interpret them as a manifestation of free topological defects. Our 
estimates of the autocorrelation exponent and the fluctuation-dissipation ratio 
differ from the XY values, while ??(Tkt) lies at the boundary of the error bar. 
Indications for logarithmic corrections at the second-order critical temperature 
are presented. However, the coupling between angles and chiralities is still strong 
and explains why autocorrelation exponent and fluctuation-dissipation ratio are 
far from the Ising values and seems stable. 



PACS numbers: 05.10.Ln, 05.50.+q, 05.70.Jk, 05.70.Ln 



Introduction 



The 2D fully-frustrated XY model (FFXY) has attracted the attention of several 
communities. Indeed, antiferromagnets on a triangular lattice with a planar anisotropy 
enter in the class of FFXY models. The FFXY also describes the superconducting- 
to-normal phase transition in Josephson-junction arrays in a perpendicular magnetic 
field for half a quantum of flux through each plaquette lj. In this work, we consider 
the following Hamiltonian 




where the sum extends over nearest neigbouring sites on the square lattice and (xi, yi) 
are the lattice coordinates of the i-th site. This Hamiltonian corresponds to the so- 
called Pileup-Domino ordering of the ferromagnetic and anti-ferromagnetic couplings. 
Frustration arises from the existence of three ferromagnetic and one anti-ferromagnetic 
couplings in each plaquette. 
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Figure 1. Example of ground state of the FFXY. The thin lines are ferromagnetic 
bonds while the bold ones are antiferromagnetic. The plus and minus signs at the 
centers of the plaquettes are the chiralities. 



The phase diagram and critical properties of the FFXY has been subject of debate 
for the last three decades. It is now established that the system undergoes two phase 
transitions: a Kosterlitz-Thouless topological transition [5] associated to the breaking 
of the invariance under the global rotation of all spins {U(l) symmetry group) and a 
second-order phase transition between a paramagnetic phase and a low-temperature 
phase where chiralities of the spins around each plaquette are anti-ferromagnetically 
ordered. This second transition is associated to the breaking of a Zi symmetry (re- 
versal of all chiralities) like in the Ising model. The existence of two phase transitions 
has been supported by many Monte Carlo simulations since the pioneering work of 
Teitel and Jayaprakash [3]. The same conclusion has been drawn for the J\ — J2 XY 
model [4]. Theoretical arguments also predict a disentanglement of U(l) and Z2 de- 
grees of freedom. Using a Hubbard-Stratonovich transformation, Choi and Doniach 
showed the equivalence of the FFXY with two coupled 0(2)-models [5]. The coupling 
is relevant under Migdal-Kadanoff renormalization and in the strong coupling limit, 
the model reduces to a Ising- XY model (the Ising spin being related to the relative 
orientation of the two 0(2)-models). Such a model have also been obtained by a 
mapping of the FFXY onto a 19- vertex model [6]. A modified version of the FFXY, 
introduced by Villain, has allowed for the clarification of the nature of the topological 
defects causing the Kosterlitz-Thouless transition. By integrating out the spin-waves 
to leave only the contribution of the topological defects, this model reduces to a gas 
of charged particles (q — ±1/2) in a uniform medium of charge 1/2 interacting via a 
logarithmic Coulomb potential [8] . 

A point of debate in the literature is whether the two phase transitions occur at 
the same temperature or at two very close ones. The breaking of the ^-symmetry 
is associated to a second-order phase transition [9] caused by the proliferation of do- 
main walls between domains of different antiferromagnetic orderings of the chirality. 
As predicted by Korshunov [10] and then observed numerically [11], at the kinks of 
these domain walls lies a q = ±1/2 topological defect and the formation of dipoles 
of two opposite charges q = +1/2 and q = —1/2 stabilizes the domain wall. As a 
consequence, the proliferation of these walls occur at higher temperatures. Recent 
Monte Carlo simulations confirm that the two phase transitions occur at two distinct 
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temperatures (see |12) and references therein). 

Another point of debate concerns the universality class of the second-order phase 
transition. Since the broken symmetry is Z(2), this universality class may be the same 
as the 2D Ising model. This idea is supported by the fact that the coupling in the 
Ising-XY model has been shown to be irrelevant at the critical point [6]. However, 
most of the Monte Carlo simulations for the FFXY or the Ising-XY model conclude 
that v ~ 0.8 while v = 1 for the Ising model (for a review of Monte Carlo estimations 
of the critical exponents, see Table 1 of ref. [H])- The most recent Monte Carlo 
simulations [T3] HU fTS] [12] now show evidences in favor of the Ising universality class. 
The disagreement with older simulations may be due to the fact that the exponents 
tend to the expected Ising values only for large lattice sizes (TSJ [12] . 

In this paper, we are interested in the out-of-equilibrium behavior of the FFXY. 
The short-time regime has been extensively studied [TH [T7] [TS] HH] in the aim to 
determine the static critical exponents as well as the dynamical exponents z and 
8. The growth of the correlation length, both for the angles and the chiralities, 
has also been studied by numerical integration of the Langevin equation [20] |2"T] . 
Temperature-dependent dynamical exponent were obtained for both angles and 
chiralities, in contradistinction to short-time dynamics simulations and relaxation of 
one-time quantities [22] for which z remains close to z ~ 2. In the following, we 
present numerical evidences that both angles and chiralities display simple aging like 
homogeneous ferromagnets. This implies in particular that the dynamical exponent 
for the angles is z = 2 for any temperature, like in the XY model. As usual in 
the context of aging, we study the two-time correlation and response functions from 
which we define a fluctuation-dissipation ratio measuring the degree of violation of 
the fluctuation-dissipation theorem. We consider two different protocols: in the first 
section, the system is prepared in one of its ground state and then let evolved below 
or at the Kosterlitz-Thouless temperature. The scaling of the two-time functions are 
compared to the XY model. In the second section, the system is prepared in the 
paramagnetic phase and then quenched at the Kosterlitz-Thouless temperature. In 
the third section, the aging of the chirality is studied when the system is prepared 
initially in the paramagnetic phase and then quenched at the Zi critical temperature. 

1. Spin- wave excitations at low temperatures 

1.1. Spin- spin correlation function in the low-temperature phase 

The system is initially prepared in one of its ground states (figure [T]) and then let 
evolved in the critical phase, i.e. at T < Tkt — 0.4461(2) [12] . with the Glauber 
dynamics [23 . If the coupling with the chirality can be neglected out-of-equilibrium, 
the thermal excitations associated to the angles are spin waves that should develop in 
the system in the very same way than in the XY model but with an effective coupling 
constant J/ \/2 [5] . Previous Monte Carlo estimates of the dynamical exponents found 
to be close to z = 2 [19] (denoted Z\ in this paper) is compatible with this picture. 
The two-time spin-spin correlation function of the FFXY are thus expected to decay 
as in the XY model [24] 
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in the spin-wave approximation. The first term is due to short-time reversible 
processes for which invariance under time translation is manifest in the depen- 
dence on t — s only (quasi-equilibrium regime). The second term is the con- 
tribution of long-time irreversible processes causing aging. It depends only on 
t/s. We have computed numerically these correlations for different final temper- 
atures T — 0.2,0.4,0.6 and 0.8Tkt- We have considered the waiting times s = 
200, 400, 600, 800, 1200, 1600, 2400, 3200, 4800, 5600, and 6400. The time t runs from s 
to 60.000. The lattice size is 256 x 256 and the data have been averaged over 12,000 
histories. 

In a first step, we have interpolated the quasi-equilibrium regime (t — s) -1 '/ 2 to 
extract the static exponent n. To eliminate the aging part, we analyzed the data along 
the curves 

- "t_f = k = Cste (3) 

VTs 

On these curves, the correlation function C(t,s) decay as (t — s)~ v / 2 so that an 
estimation of n is easily obtained by a linear fit of In C (t, s) versus \n(t — s) . Practically, 
this was accomplished by first doing an interpolation of the correlation function in the 
plane (t,s) between the simulation points. On figure EJ our estimates of n/2 are 
plotted versus l/n. Since we are considering the quasi-equilibrium regime, t ~ s, 
the value of ij has to be measured in the regime of small k and thus on the right of 
the figure. Deviations, visible on the left part of the figure, are due to corrections to 
the dominant behavior of the correlation function. To test for the possible existence 
of transient regime at small s, we have also tried to remove the smallest waiting 
times: the different symbols on the figure correspond to effective exponents taking 
into account all waiting times (circles), only waiting times s > 400 (squares), s > 600 
(diamonds) up to s > 5600 (stars). Even though the statistical errors increase when 
fewer waiting times s are taken into account, all curves are nicely compatible. Our final 
estimates of r\ are given in Table [T] They are in good agreement with other estimates 
found in the literature (see figure [3]), including static Monte Carlo simulations, and 

close at low temperature to the linear behavior r\ — k °^P' predicted in the spin-wave 
approximation. 



T/T c 


rj/2 (MC) 


77/2 (SW) 


0.2 


0.0106(6) 


0.0100 


0.4 


0.0225(9) 


0.0201 


0.6 


0.0358(8) 


0.0301 


0.8 


0.0536(20) 


0.0402 



Table 1. Comparison of the Monte Carlo estimates (MC) of the exponent r)/2 
obtained in the quasi-equilibrium regime with the predicted values in the spin- 
wave approximation (SW). 



The exponent 77 can be estimated from the aging part of the correlation function 
C(t,s) too. Since it is necessary to take into account correlation functions at well 
separated times s and t, the correlation function is very small and thus noisy. A 
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Figure 2. Effective exponent r//2 obtained by interpolation of the spin- 
spin correlation function in the quasi-equilibrium regime. The different curves 
correspond to temperatures 0.2Tkt, OATkt, O-GTkt and 0.87Vt (from bottom 
to top). The different symbols correspond to different interpolation ranges from 
s > 200 (circle) to s > 5600 (star). 
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Figure 3. Comparison of the exponent r\ with other estimates found in the 
literature and the prediction of the spin-wave approximation. 
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Figure 4. Effective exponent ?//2 obtained by interpolation of the spin-spin 
correlation function in the aging regime. The different curves correspond to 
temperatures 0.2T C , 0.4T C , 0.6T C and 0.8T C (from bottom to top) and the different 
symbols to different waiting times s. The exponents r;/2 obtained in the quasi- 
equilibrium regime are plotted as a dashed line. On the right, the same data are 
plotted with respect to l/t m i n . 



precise determination of n is more difficult than in the quasi-equilibrium regime. 
As a consequence, we will restrict ourselves to check that the previously estimated 
exponents n are compatible with our data in the aging regime. To that purpose, we 
made a log-log interpolation of (t — s) I? / 2 C(i, s) versus ^j=. To take into account 
possible corrections to the asymptotic behavior, we varied the interpolation window. 
Effective exponents are plotted on figure|4]with respect to the smallest value of x = ^fc= 
entering into the fit. The rightmost points are obtained by an interpolation over all 
times t 6 [s + 1; 60, 000]. The leftmost points come from an interpolation over only 
the three last times t £ [59, 998; 60, 000]. The different symbols correspond to different 
waiting times s. The effective exponents display a plateau for the values obtained in 
the quasi-equilibrium regime and plotted as a dashed line on figure [4] However, as 
temperature increases, the width of this plateau shrinks, ft could tempting to interpret 
the deviations from the plateau as the effect of corrections to the aging behavior x v ^ 2 . 
However, as can be seen on the right figure [H where the effective exponents are now 
plotted with respect to the smallest time t entering into the fit, these deviations depend 
only on t and occur for large values of t. They are thus due to an under-sampling 
of the correlation at large times t that induces large systematic fluctuations of the 
effective exponents when only the last few points enter into the interpolation. They 
are thus numerical and not physical. 

In the case of a quench at the Kosterlitz-Thouless temperature Tkt, the 
estimation of the exponent n using the interpolation of C(t,s) ~ (t — s)~ ?J / 2 in 
the quasi-equilibrium regime becomes much less accurate. We obtain an estimate 
7/ ~ 0.27(2). However, a nice collapse of the scaling function (t ~ s) v ^ 2 C(t : s) versus 
the scaling variable (t + s)/\/~ts is observed as shown on figure [5] and gives evidences of 
the validity of the expression ([2]) up to the Kosterlitz-Thouless transition temperature. 
The effective exponents do not present any plateau but only seem to tend towards the 
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Figure 5. Scaling function (t — s) v ^ 2 C(t, s) of the spin-spin correlation function 
versus the scaling variable (t+s) /vis for the FFXY model prepared in the ground 
state and quenched at its critical temperature. The different curves correspond to 
different waiting times (s = 200, 400, 600, 800, 1200, 1600, 2400, 3200, 4800, 5600, 
and 6400). The dashed curve is the expected asymptotic behavior \{t + 

,)/v^]" /2 . 



expected value. 



1.2. Response function in the low-temperature phase 

The linear response of the observable to, say the magnetization, at time t to an external 
field h, say the magnetic field, coupled to the system at the instant s is defined as 

».)=(«) (4) 

At equilibrium, this quantity is related to the correlation function C(t, s) = 
(m(t)m(s)) — (m(t))(m(s)) by the fluctuation-dissipation theorem (FDT) 

k B TR(t,s) = -^-C(t,s) (5) 
os 

The violation of this theorem in aging systems has attracted much attention these last 
ten years. In the context of mean-field spin-glasses, the following extended relation 
was introduced [25] 

k B TR(t,s)=X(t,s)-^C(t,s) (6) 

where X(t,s) differs from the value X = 1 out of equilibrium. In the aging regime, 
X(t, s) is believed to tend towards a universal quantity, denoted I^,. 
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In the case of the FFXY model, the local order parameter is the two-dimensional 
vector As a consequence, correlation and response functions are tensors: 

C^(t,s) = ((^(t).^)(ffi(«).^)) - ((a i (t).e M ))((o ! i (s).e;)) (7) 

R ^ s) = lh^sJ- e » (8) 

where the magnetic field is decomposed ash — h v e u . In the following, we will consider 
only the traces Tr C(t,s) and Tr R(t,s) corresponding to the usual expression of the 
correlation function: 

Tr C(t,s) = (cos0;(i) cos0;(s)) + (sin0;(t) sin^(s)) 

= (cos(^(i) - 0i{8))) (9) 

where s,; = ( cos 9i sin 9i ) . This choice was already made in the study of the 
fluctuation-dissipation ratio in the XY model [26j . 

Since the FFXY behaves at low temperature as an XY model with a stiffness 
J/y/2, the response function can be calculated in the spin-wave approximation. This 
leads to [57] 

RQt,s)~{t-8)- 1 -*' 2 ( t -±Z) (10) 



ts . 

The response function can be computed directly in Monte Carlo simulations [29] . 
In principle, it could be used to estimate n from (fTU]) . In practise, the response 
function is noisier than the correlation function and thus does not allow for a precise 
determination of rj. However, with a sufficiently large number of histories, the 
fluctuation-dissipation ratio X(t,s) can be computed with a reasonable accuracy. 
Using the scaling behavior for the correlation function ([2]) and the response function 
(|10| in the spin-wave approximation, one can calculate the scaling behavior of the 
fluctuation-dissipation ratio: 

v . k B TR(t,s) s(t + s) 

X(t,s)- 



d s C{t, s) t 2 - Us - s 2 

that displays a divergence for t/s = 2+y5. Such a divergence was already reported [28] 
for the XY model. We also observe a divergence of X(t,s) for the FFXY but its 
location tends to increase with the waiting time s. Unfortunately, the computation of 
the fluctuation-dissipation ratio is very CPU time-demanding and we had to restrict 
ourselves to a lattice size 192 x 192 and a final time tt = 4000. The data have been 
averaged over 15,000 histories at the two lowest temperatures and 45,000 for the 
others. The waiting times are s — 60, 80, 100, 120, 140, and 160. As shown on figure 
[7] the location of the divergence io/so displays a nice linear behavior with 1/s with a 
value at 1/s — * compatible with the predicted value 2 + y5 ~ 4.236. This provides 
an indirect confirmation of the scaling behavior (flQ]) . 
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Figure 6. Fluctuation-dissipation ratio X(t,s) at T = 0.2Tkt- The different 
curves correspond to different waiting times s. 
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Figure 7. Location to /so of the divergence of the fluctuation-dissipation ratio 
X(t, s) versus l/s. The different curves correspond to different temperature. The 
dashed line is the expected value in the spin-wave approximation. 
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2. Aging of angles 

We now consider the following protocol: the system is prepared in a random state, 
i.e. at infinite temperature, and then quenched at the Kosterlitz-Thouless transition 
temperature. As we will show, the dynamics of the FFXY shows aging in a way 
similar to the XY model. Above the Kosterlitz-Thouless transition temperature, the 
spin-wave approximation is not expected to hold anymore but we can assume, like in 
the XY model, that the two-time spin-spin correlation function C{t, s) decay as 

C(M)~(«-.)-^/(f)~(.-.)^/c(^f) (U) 

The prefactor (t— s)~ v / 2 corresponds to the quasi-equilibrium regime while the scaling 
function / encodes the aging behavior of the correlation function. According to the 
aging hypothesis, the relevant scaling variable is £(i)/£(s) where the characteristic 
length of domains is expected to grow as £(i) ~ (t/lnt) 1 /*, i.e. algebraically 
(infinite relaxation time) with a logarithmic correction due to the pinning of the 
domain walls by the free vortices present in the paramagnetic phase [30) . In the 
case of the XY model, the dynamical exponent is z = 2. To test eq. (TTTT) . we 
have performed Monte Carlo simulations for a lattice size 256 x 256 with a final 
time tf = 60,000. The data have been averaged over 5000 histories. The data are 
noisier than in the low-temperature and it is not possible to obtain a more accurate 
estimate of r\ from the quasi-equilibrium regime. On figure [SJ the scaling function 
fc = (t— s) v ' 2 C(t, s) is plotted versus ^4- The different curves correspond to waiting 
times s = 500,1000,2000,4000,8000, and 12000. The collapse of the different curves 
has been obtained using a c — 0.1. We are left with an algebraic decay fc(%) ~ that 
confirms the scaling hypothesis (|11|). A close inspection of the effective exponents for 
the different waiting times reveals an evolution with s. On the inset of fig[8]the effective 
exponents calculated using a power-law fit fc{x) ~ in the window [t m i n ;i/] are 
plotted versus lni m i„/i m i n . The different curves correspond to the different waiting 
times s. Again, the large fluctuations for the leftmost points (larger values of t) occur 
at the same value of t, independently of the waiting time s. We thus interpret them 
as due to an under-sampling of the correlation causing systematic deviations of the 
effective exponents when the number of points entering into the interpolation becomes 
small. As a consequence, we have to consider the plateaus. An evolution with s of the 
effective exponents on these plateaus is clearly perceptible. The exponents are plotted 
on figure [8] The extrapolation for s — >• +00 leads to <j> ~ 0.74 and thus, defining the 
autocorrelation exponent as [32] C(t,s) ~ t~ x l z ~ £~' ? / 2 ~^, we obtain X/z ~ 0.84. 
This value is slightly larger than the one obtained in the short-time regime from the 
initial-slip exponent 6 (X/z = | — 6 = 0.808(3)) [15]. Note that short-time dynamics 
is equivalent to our analysis with s = [33) and is thus unable to give a correct 
estimate if the effective exponent evolves strongly with the waiting time s. Both our 
estimate of X/z and that obtained by Short-Time Dynamics are incompatible with 
the autocorrelation exponent of the XY model: X/z = 0.625 [19], 0.738(4) [22] or 
0.738(5) [34]. 

Again the response is too noisy to investigate its asymptotic scaling behavior but, 
using smaller lattice sizes and final times, it is still possible to improve significantly 
the statistics and estimate accurately the fluctuation-dissipation ratio. We have used 
a lattice size 192 x 192 with time running up to tf = 4, 000. This allowed us to average 
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Figure 8. On the left, scaling function (t - s)' ) '' 2 C(t, s) versus of the 

spin-spin correlation function for the FFXY prepared in a random state and 
quenched at the critical temperature. The different curves correspond to different 
waiting times (s = 500, 1000, 2000, 4000, 8000, and 12000). The inset corresponds 
to effective exponents calculated using a power- law fit in the window [i m i n ; tr] and 
plotted versus lnt m ; n /t m ; n . On the right, the effective exponents on the plateau 
are plotted versus 1/ ' yfs. The red line is a linear fit of the data. 



the data over 500, 000 histories. On figure El the FDR is plotted with respect to jj^f ■ 
A good collapse of the data for different waiting times is observed. This collapse is 
not as good if logarithmic corrections are omitted. The asymptotic value of the FDR 
is estimated to be X^ — 0.385(15). In contradistinction to the algebraic decay of the 
spin-spin correlation function, X^ does not present any evolution with the waiting 
time s. For comparison, we have performed a similar simulation for the pure XY 
model. The asymptotic FDR is estimated to be = 0.215(15) for the XY model, 
i.e. very different from our value for the FFXY. 



3. Aging of the chirality 

We have also investigated the dynamics of the chirality, denoted \ m the following. 
For each plaquette, the local chirality is defined as 

Xn = sign [ ^ Jy sin(0j - 0,) (12) 

where the sum extends over the four bonds forming the plaquette. Chirality is 
expected to behave as antiferromagnetically-coupled Ising spins but a controversy 
exists regarding the belonging of the associated transition to the Ising universality 
class. In the following, we will compare two-time chirality-chirality correlation 
functions and the FDR associated to the response to a field coupled to the chirality 
with those of the Ising model. 
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Figure 9. Fluctuation-dissipation ratio associated to the angles, for the FFXY 
prepared in a random state and quenched at the critical temperature. The 
different curves correspond to different waiting times (60, 80, 100, 120, 140, 160, 
and 180). In the insert, the FDR is plotted for the XY model under the same 
protocol. The different curves correspond to s = 10, 20, 50, 100, 200, and 400. 




Figure 10. Chirality-chirality correlation function C{t, s) versus tlns/slnt for 
the FFXY prepared in a random state and then quenched at the inverse critical 
temperature T c = 0.4545. The different curves correspond to different waiting 
times (500, 1000, 2000, 4000, 8000, and 12000). In the insert, the effective exponent 
X/z computed by a power-law interpolation (t — s) ac C(t, s) ~ (t/ \nt)~ x / z+a " 
between t and tf is plotted with respect to [tlns/slnt] . The dashed line 
corresponds to the exponent X/z for the pure model. 
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The system is initially prepared in the paramagnetic phase and then quench at 
the Z 2 critical temperature. We have considered two different estimates of the critical 
temperature: T c = 0.4545(2) [TTJ HH] and T c = 0.453243(2) [TO]. We have computed 
the two-time chirality-chirality correlation function 



for a lattice size 256 x 256 up to the final time tf — 200, 000. Data have been averaged 
over 5000 histories. During a quench at the critical temperature, autocorrelation 
functions are expected to behave as C(i, s) ~ (t — s)~ Qc /(£(i)/£(s)) where the first 
term corresponds to the quasi-equilibrium regime characterized by a critical exponent 
a c = 2(3/ vz c and the second one is the aging part depending only on the ratio £(i)/£(s). 
In the FFXY, we assume that the existence of free topological defects in the high- 
temperature induces logarithmic corrections in the growth law. Thus the correlation 
length for the chirality-chirality correlations is expected to grow as £(i) ~ (t/ lni) 1 ' 20 , 
i.e. in the same way as for the spin-spin correlation functions. As a consequence, we 
expect the scaling behavior 



This behavior of the chirality-chirality correlation function is observed numerically 
only with the estimate T c — 0.4545 of the critical temperature. Figure [TO] shows the 
collapse of the scaling function fc = (t — s) ac C{t, s) versus t\ns/slnt for different 
waiting times. This collapse is clearly destroyed when a c is chosen smaller than 0.08 
or larger than 0.12. This leads us to a c = 0.10(2), compatible with the Ising value 
a c = 2(3 /vz c ~ 0.115. We note that removing the logarithmic corrections does not 
affect significantly the quality of the collapse. For the second and most accurate esti- 
mate of the critical temperature T c = 0.453243(2), no collapse of the scaling function 
is obtained for any positive value of a c (only a slightly negative value would bring the 
data to collapse). This shows that the system is not critical in the quasi-equilibrium 
regime but already in the low-temperature phase. In the following, we will consider 
only the estimate T c = 0.4545 of the critical temperature. 

As shown on figure [TO] the aging part fc(x) of the chirality-chirality correlation 
function decays algebraically. To recover the usual form C(t,s) ~ s~ ac (t/s) X ^ z in 
the asymptotic regime, the scaling function fc{x) needs to decay as x ~ x / z+ac . In 
the inset of figure [TO] is plotted the effective exponent X/z computed by a power-law 
interpolation in the interval [t,tf] where t runs from s to tf — 2. This effective ex- 
ponent is plotted with respect to fp^ = x for the different waiting times s. The 
exponent X/z is recovered in the limit x — * 0, i.e. at the left side of the graph. Of 
course, in this region the number of points entering in the interpolation is small and 
thus the accuracy of the effective exponent decreases. Moreover, the data have been 
produced by a slow dynamics and are thus correlated. The smooth variations at the 
left are really fluctuations and not any physical effect. From this plot, we estimate 
X/z — 0.98(5), i.e. incompatible with the Ising value X/z ~ 0.738 (dashed line on 
the graph). Even though the curve of the effective exponent for the largest waiting 
time is above the others, this cannot be interpreted as an evolution of the exponents 
with s. Indeed, the sequence of curves is not ordered according to the waiting times. 
Our conclusion is remarkably different from that drawn for the Ising-XY model in the 
low-temperature phase for which an estimate A ~ 1.25 compatible with that of the 



C(t,s) = ( X (t) X (s)) 



(13) 




(14) 
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pure Ising model was measured already for a relatively small final time tj = 1000 [31J. 
For completeness, we should note that when we remove the logarithmic corrections, 
the effective exponent is smaller, X/z ~ 0.90(5) but still incompatible with the Ising 
value. However, we believe that logarithmic corrections are necessary as we will see 
in the following when considering the scaling of the fluctuation-dissipation ratio. 

The fluctuation-dissipation theorem is also violated at the Z2 phase transition. 
A Zeeman Hamiltonian 



□ 

that couples the chirality of each plaquette to an external field ha, allows for the 
definition of a response function 



This quantity is evaluated numerically |29J . We can then compute the fluctuation- 
dissipation ratio X(t, s) = ksTR(t, s)/d s C(t, s). We used a lattice size 192 x 192 with 
a final time tf = 2000. Data have been averaged over 10 6 histories. The FDR is 
plotted on figure [TT1 with respect to s/t (left) and slni/tlns (right). A better collapse 
is obtained with logarithmic corrections, showing that the dynamics of the domain 
walls between different antiferromagnetic orderings of the chirality is sensible to the 
vortices. The asymptotic value Xoo is estimated to X^ = 0.405(5) (a compatible 
value would be obtained if logarithmic corrections would not be taken into account). 
This value is significantly different from the Ising value Xoo ~ 0.328 [35]. Again, no 
evolution of 1^ with the waiting time s is observed. 

Conclusions 

We have given numerical evidences that the thermal excitations in the low-temperature 
phase are spin- waves whose dynamics is the same as in the XY model. The dynamical 
exponent is thus z = 2 in the whole critical phase. The exponent n that can be 
estimated from the spin-spin two-time autocorrelation function is compatible with 
the prediction of the spin-wave approximation and with estimates from static Monte 
Carlo simulations close to Tkt- The fluctuation-dissipation ratio displays divergences, 
like in the XY model, at some locations that tend toward the value predicted by the 
scaling theory. When initially prepared in the paramagnetic phase and then quenched 
at Tkt (resp. T c ), spins (resp. chiralities) age in the same way as an unfrustrated 
ferromagnet. Spin-spin autocorrelation function decays slightly faster (X/z ~ 0.84) 
than in the XY model. The FDR takes asymptotically a very different value too. 
These results are consistent with the fact that our estimate of n = 0.27(2) is greater 
than the XY value, even though at the boundary of the error bar. The chirality- 
chirality autocorrelation function decays algebraically at the critical point with an 
exponent X/z ~ 0.98(5) incompatible with that of the Ising model. The evolution 
of the effective exponent with the waiting time was taken into account. The FDR 
is also fully incompatible with the Ising value. No evolution of these values with 
the waiting time is visible, in contradistinction to static Monte Carlo simulations for 
which a cross-over is clearly observed when increasing the lattice size [T3J [T2] . The 
coupling between angles and chiralities is probably still too strong in the range of 




(15) 




(16) 
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Figure 11. Fluctuation-Dissipation Ratio X(t,s) associated to the chirality 
versus t/s (left) and tlns/shit for the FFXY prepared in a random state and 
quenched at the critical temperature. The different curves correspond to different 
waiting times (fO, 20, 30, 40, 50, 60, and 80). The dashed line corresponds to the 
asymptotic value for the Ising model. 

parameters used in our simulations. This coupling is manifest in the fact that the 
FDR is affected by logarithmic corrections that may be interpreted as an influence 
of the topological defects causing the Kosterlitz-Thouless transition. Moreover, the 
fact that our estimate of X^ for the chirality is relatively close to the value measured 
at the Kosterlitz-Thouless temperature for the angles may also be the result of a 
strong coupling. To reach the cross-over regime and eventually enter the true critical 
regime, the correlation length should be allowed to reach much larger values. This 
requires of course larger lattice sizes but also larger times since £ ~ (t/ mi) 1 / 2 . The 
computational effort scales with the lattice size as L 2+z . It is therefore identical to 
static Monte simulations since the number of Monte Carlo steps needs to be increased 
when the lattice size increases and only local updates are possible for the FFXY. 
However, it seems that the amplitude is playing against us! 
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